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Abstract 

MCViNE (Monte-Carlo Virtual Neutron Experiment) is an open-source Monte 
Carlo (MC) neutron ray-tracing software for performing computer modeling and 
simulations that mirror real neutron scattering experiments. We exploited the 
close similarity between how instrument components are designed and oper¬ 
ated and how such components can be modeled in software. For example we 
used object oriented programming concepts for representing neutron scatterers 
and detector systems, and recursive algorithms for implementing multiple scat¬ 
tering. Combining these features together in MCViNE allows one to handle 
sophisticated neutron scattering problems in modern instruments, including, 
for example, neutron detection by complex detector systems, and single and 
multiple scattering events in a variety of samples and sample environments. In 
addition, MCViNE can use simulation components from linear-chain-based MC 
ray tracing packages which facilitates porting instrument models from those 
codes. Furthermore it allows for components written solely in Python, which 
expedites prototyping of new components. These developments have enabled 
detailed simulations of neutron scattering experiments, with non-trivial sam¬ 
ples, for time-of-flight inelastic instruments at the Spallation Neutron Source. 
Examples of such simulations for powder and single-crystal samples with various 
scattering kernels, including kernels for phonon and magnon scattering, are pre¬ 
sented. With simulations that closely reproduce experimental results, scattering 
mechanisms can be turned on and off to determine how they contribute to the 
measured scattering intensities, improving our understanding of the underlying 
physics. 
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spectrometry 


1. Introduction 

Data analysis for neutron time-of-flight spectroscopy has predominantly been 
a conversion of raw data to intensities with instrument independent units; 
namely S(Q,ui) or S(Q,E), where S is the dynamical structure factor, and 
Q, w, E are the wave vector, frequency, and energy. This approach has been 
the basis for many studies, but some cases require a more elaborate analysis 
scheme (see, for example, [1, 2, 3, 4, 5, 6, 7, 8, 9]). 

More specifically, traditional reduction assumes single scattering, that the 
resolution function can be modeled by a simple function, that the scattering 
process of interest can be readily separated from other scattering processes in 
the system, and that scattering from the sample environment is negligible or 
easily subtracted. Monte Carlo ray-tracing can provide a straightforward way 
to perform analysis in cases where some, or all, of the above assumptions are 
not valid. Specifically it can be used to investigate complex samples whether 
they have multiple scattering, multiple modes of scattering, i.e. magnetic and 
lattice vibrations, or consist of a conglomeration of scatterers like a sample and 
sample environment. Additionally as the instrument is modeled component by 
component, instrument resolution is inherently included in the simulation. 

As the promise of Monte Carlo ray tracing was clear, it has been used in 
codes that handle a subset of these complexities. MCS[10] and MSCAT[11] were 
early packages that used Monte Carlo methods to compute multiple scattering. 
Later, several general-purpose Monte Carlo neutron instrument simulation pack¬ 
ages were developed and optimized to help design neutron instruments, namely 
McStas [12, 13], VITESS [14], IDEAS [15], and NISP [16]. Less-general-purpose 
MC programs exist in some popular neutron data analysis software packages, 
including DISCUS [17, 18, 19] and RESTRAX [20, 21], In comparison, relatively 
few studies [22, 23, 24, 25, 26] have employed full-fledged Monte Carlo ray trac¬ 
ing to help analyze experimental results from neutron scattering measurements, 
while there are growing efforts to perform virtual neutron experiments using 
MC ray tracing [27, 28, 29, 30, 31]. MCViNE[32], an open source software, is 
designed to easily allow complex studies (see [9] for an example) and therefore 
should accelerate the use of Monte Carlo ray tracing in the analysis of experi¬ 
mental data. 

The improvements in computing over the last decades include hardware de¬ 
velopments, the emergence of object-oriented languages C++ and Python, and 
advances in software design [33]. The amazing increases in hardware capabilities 
now allows detailed simulations to run in a reasonable time. The impact of mod¬ 
ern software engineering practice on MC simulation codes is still progressing. 
MCViNE was developed, as part of the DANSE project [34], using such soft¬ 
ware practices. Therefore it uses object oriented programing (OOP) concepts 
to represent instrument and sample components, following the hierarchical and 
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modular nature of the neutron instruments, detector systems, and sample en¬ 
vironments. This use of hierarchies from inside an instrument component or 
sample kernel up through the total instrument goes beyond the imperative pro¬ 
gramming paradigm, which is popular in other codes. 1 By taking advantage 
of OOP design patterns as well as recursive algorithms, the MCViNE archi¬ 
tecture supports an easily extensible library for scattering kernels suitable for 
both samples and detector systems. This approach allows for maximum flexibil¬ 
ity, extensibility, and reuse of scatterer arrangements, geometrical shapes, and 
scattering mechanisms, and hence improves the sustainability of the software. 

This paper contains a description of the MCViNE software framework. Sec¬ 
tion 2 presents the challenges in simulating experiments carried out in modern 
neutron instrumentation, and an overview of the software engineering design 
and main software components and algorithms. Examples of MCViNE simu¬ 
lations for powder and single-crystal experiments are presented in Section 3. 
Concluding remarks follow in Section 4. 

2. Simulation of Spectra Measured at Modern Neutron Spectrome¬ 
ters with MCViNE 

2.1. Challenges for simulating neutron scattering experiments 

Monte Carlo (MC) ray tracing simulations of neutron scattering spectrom¬ 
eters with support of multiple scattering were performed from 1970s using 
MCS[10] and MSCAT[11]. They were used to understand the effects of mul¬ 
tiple scattering on the measured spectra [35]. In the 1990s, with increasing 
needs of simulating neutron instruments for the purpose of instrument design 
and optimization, several MC neutron ray tracing packages emerged, including 
McStas[12, 13], Vitess[14], Ideas[15], and NISP [16]. Simulations performed with 
these packages provide not only independent checks for analytical calculations 
of instrument performance, but also can be used to obtain optimal parameters 
for instrument design (see, for example, [36, 37, 38, 39, 40, 41, 42, 43, 44, 45]). 
Unlike MCS and MSCAT, most of these newer MC software packages (with a 
notable exception in NISP) treat simulation of a neutron instrument as a linear 
chain of neutron optical components, each of which modify some set of neutron 
beam characteristics such as spatial divergence and energy distribution. This 
linear approach greatly improved the computational efficiency and simplified 
the coding of instrument simulations. Such an approach is adequate because 
the physical formation of a neutron scattering instrument is linear, and the 
underlying hypothesis that neutrons at the downstream components are rarely 
scattered back to the upstream components holds up well in most instrument 
configurations. As a result, these software packages, especially McStas and 
Vitess, are making significant impacts for neutron instrument design. 

For Monte Carlo neutron ray tracing simulations to be useful for interpre¬ 
tation of neutron scattering spectra, it is necessary to include detector systems 


1 A brief discussion is available in the supplemental material. 
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and samples/sample environments, for which the physical arrangements are of¬ 
ten non-linear. The detector systems in modern neutron scattering instrumen¬ 
tation display modularity, repetition, and sometimes hierarchical organization, 
reducing the engineering difficulties in manufacturing, testing, and validation. 
For example, the four direct geometry time-of-flight spectrometers at Spalla¬ 
tion Neutron Source (SNS) [46], ARCS [47], SEQUOIA [48, 49], CNCS[50], and 
HYSPEC[51] share some instrument elements[52]: the Fermi chopper slit pack¬ 
ages for the ARCS and SEQUOIA instruments are interchangeable, and the 
detector systems of all four instruments use the same so-called 8-packs [53], 
each of which is a detector pack consisting of eight 3 He linear position sensi¬ 
tive detector tubes (LPSD). The 8-packs are arranged in a vertically-oriented 
cylindrical geometry around the sample position, forming a hierarchical orga¬ 
nization of pixels, tubes, packs, and detector rows. A simplified illustration of 
the detector hierarchy can be found in Figure 1(c). 

Samples, sample holders, and sample environments are constituents of a 
“sample assembly”, the term we use to refer to the collection of neutron scat- 
terers near the sample position. A sample-assembly example is illustrated in 
Figure 1(b). A linear representation faces challenges in simulations of a sam¬ 
ple assembly because it does not match physical arrangements of a sample or 
samples with associated sample environment components, and because neutrons 
can often scatter back and forth among the constituents. A good example of 
complex sample environments is the MICAS furnace [54]: the incident neutron 
beam intersects heating elements and up to 8 heat shields, all thin vertical hol¬ 
low cylinders, in addition to the sample. The sample itself can be challenging 
too. For example, “single crystal” superconductor samples for inelastic neu¬ 
tron scattering experiments often consist of co-aligned arrays of small crystals 
mounted on an aluminum plate [55]. Another major difficulty stems from the 
fact that a variety of different scattering mechanisms may be present in the 
combination of sample, sample holder, and sample environment. 

2.2. General capabilities of MCViNE 

MCViNE[32] is a general purpose neutron ray tracing package that combines 
the two different approaches taken by MSCAT-like packages and McStas-like 
packages. At the instrument level, it allows users to create a simulation appli¬ 
cation as a linear chain of neutron components, each fully configurable by its 
type and corresponding geometrical and physical properties using either com¬ 
mand line options or an xml-based configuration file. At the component level, 
two general components exist with support of hierarchical representation for the 
sample assembly and the detector system. A detector system can be specified 
with an xml file describing its hierarchy, while a sample assembly is specified by 
a collection of xml hies, one for the geometrical organization of the constituent 
scatterers, and others describing the scattering mechanisms for them. Multiple¬ 
scattering among scatterers in a sample assembly can be turned on and off by 
a command line option. 
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2.3. Architecture Overview 

MCViNE is implemented in C++ and Python. The C++ portion contains 
fundamental mathematics, the neutron ray tracing mechanism, the mechanism 
to include components from other MC ray tracing packages, and support of 
multiple scattering and composite scattering kernels. These various functions 
of the C++ code are implemented in two layers, which are discussed in the 
supplementary information. There are two aspects of the C++ usage that 
bear more discussion. The components from McStas are included in the CH—h 
layer so that they can be used as individuals rather than as a whole compiled 
instrument. Also, the combination of the object oriented features of C++ and 
the speed of the compiled code make this the appropriate place to implement the 
composite kernel and scattering kernels. Currently two such components exist, 
one for sample assemblies and one for detector systems. This unique feature of 
MCViNE, the “composite scatterer”, will be explained in detail in section 2.4. 

A Python layer on top of the two C++ layers allows construction of a com¬ 
ponent chain similar to McStas-like neutron ray tracing packages, provides the 
interface to the C++ components to include in that chain, and allows for intro¬ 
duction of components that are completely written in Python. This last feature 
makes it extremely easy to create simple neutron components, and to create 
prototypes of more sophisticated neutron components. 

2-4■ Composite Scatterer 

This section starts with basic concepts with regard to the “composite scat¬ 
terer” , followed by an introduction to the essential object-oriented software de¬ 
signs centered around it, and finishes with examples of ray-tracing algorithms 
in sample assemblies and detector systems, enabled by these designs and needed 
in simulation of neutron scattering experiments in modern spectrometers. 

2.4-1- Concepts 

In MCViNE, a “composite neutron scatterer” represents a group of physical 
objects, for example, a powder sample in an aluminum can, a single crystal 
sample surrounded by a furnace, or a detector system. An “elemental scatterer” 
is a scatterer without constituent scatterers. A “homogeneous scatterer” is one 
kind of elemental neutron scatterer, whose scattering function is homogeneous 
within its volume. The scattering properties are modeled using one “scattering 
kernel” or a combination of several “scattering kernels”, each of which represents 
one scattering mechanism, such as incoherent one-phonon nuclear scattering or 
coherent magnetic scattering. 

These concepts are illustrated in Figure 1. Figure 1(a) is an example of 
an abstract hierarchy of a composite scatterer in MCViNE. In principle, the 
hierarchy in MCViNE can be arbitrarily deep. In practice, the depth of the 
hierarchy is limited by factors such as computing resources available for the 
simulations, and compiler limitations. Figure 1(b) depicts a sample assembly 
that is a composite of two-level hierarchy, in which the bottom level consists 
of two homogeneous elemental scatterers: one aluminum can and one copper 
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Figure 1: Concepts in scattering composites: composite scatterer, homogeneous scatterer, 
shape, scattering kernels, (a) an abstract hierarchy of an abitrary scattering composite, (b) 
and (c) are concrete examples of such hierarchies, (b) top view of a sample assembly consisting 
of an aluminum can and a copper plate, (c) a detector system - how it is constructed from 
an elemental scatterer (detector tube) in a three-level hierarchy. 
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plate sample. Figure 1(c) represents a detector system consisting of ? He eight- 
packs that form roughly a cylindrical arrangement around the sample position. 
The detector system is represented in MCViNE in a three-level hierarchy: at the 
bottom level is the 3 He detector tube; at the middle level, 8 such tubes construct 
an 8-pack; at the top-level, the detector system consists of a collection of 8-packs. 
Such hierarchical representations allow MCViNE to model the physical reality 
closely. 

Object Oriented Designs 

These concepts lead to one major design decision made in the MCViNE 
project: employment of the “composite design pattern”, which is used to de¬ 
scribe the part-whole relationship and to represent a hierarchy, and allows clients 
to treat composites and their constituents in a uniform way, and the “visitor 
design pattern”, which allows separation of operations from the objects to be 
operated on, so that new operations can be added without touching these ob¬ 
jects [33]. By using these design patterns for scattering composites, we can unify 
the programming interfaces to the operations on both the composites and the 
individual elemental objects. Composite and visitor patterns are used in three 
major aspects of the MCViNE neutron scattering model: the neutron scatterers, 
the geometric shapes of scatterers, and the scattering kernels. 

Neutron scatterers. By using the composite pattern, algorithms for mul¬ 
tiple scattering can be consolidated in one implementation. Scattering from a 
composite neutron scatterer starts with a determination of which constituent 
intersects the incident neutron ray, and then delegates the scattering assessment 
to that particular constituent, which could be a composite itself that requires 
another delegation for scattering. The hierarchical representation of neutron 
scatterers and this recursive algorithm work for both samples and detector sys¬ 
tems, and can improve computing efficiency and code maintenance. 

Geometric shapes Using constructive solid geometry (CSG) (see for exam¬ 
ple [56]), composite shapes are constructed from basic shapes such as cylinders 
and blocks, and composites by using operations such as union, intersection and 
difference. Ray-tracing through shapes is therefore simplified as visitor methods 
of the primitive shapes and the binary shape operations. 

Scattering kernels Composite scattering kernels make it easy for the Monte 
Carlo algorithm to sample a total S(Q, w) consisting of both slowly-varying re¬ 
gions and regions containing sharp features, such as diffraction and coherent 
phonon scattering, by allowing users to combine different kernels such as in¬ 
coherent and coherent kernels in a kernel composite. In addition, users can 
organize scattering kernels into groups; this makes it easy to apply importance 
sampling (assign different weights to different kernels or kernel groups) in sim¬ 
ulations. 

A “scattering kernel” in MCViNE is conceptually different from sample com¬ 
ponents in linear MC neutron ray tracing packages. A scattering kernel in 
MCViNE is an abstraction of the scattering mechanisms such as diffraction, 
nuclear scattering by phonons, and magnetic scattering by spin waves. It does 
not include the sample geometry but only the scattering physics. A sample 
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component in earlier packages, on the other hand, includes both the geometry 
and the physics in one programming unit. By separating the implementation 
of “scatterer”, “shape”, and “scattering kernel”, a sample in MCViNE can con¬ 
sist of a combination of scattering kernels. Furthermore, the scattering kernel 
library in the MCViNE framework can be extended without affecting the logic 
of geometric ray tracing, which is implemented in “shape” and “scatterer”. For 
example, for isotropic scattering, scattering kernels taking a histogram form of 
S(Q,lo) can be supported, as well as phonon scattering kernels taking phonon 
energies and polarizations as inputs (examples are given in Section 3.1 and 
3.2). A scattering kernel conveniently taking the analytical form of a dispersion 
can be used (see [9] for an example) to improve convergence (effective when 
combined with other scattering kernels), and to avoid unnecessary broadenings 
resulted from approximating S(Q,lu) using a histogram. Inelastic and elastic 
scattering kernels for single crystals can also be developed to simulate single 
crystal experiments. 

2-4-3. Algorithms 

In this section, the general ray tracing procedure is first described briefly, and 
then implementation of multiple scattering and ray tracing in detector systems 
are presented as algorithm examples that benefit from the conceptual analysis 
and the software design. 

Ray Tracing 

In components such as sample assemblies or detector systems, neutron scat¬ 
tered are represented by a hierarchy of objects with shapes and scattering mech¬ 
anisms. Ray tracing of a neutron happens by first investigating which one of the 
neutron scatterers at the top level of the scatterer composite hierarchy inter¬ 
cepts the neutron. This is done by computing the intersections of the forward 
ray of the neutron and all the shapes of the top level constituents. A random 
selection might be necessary if multiple top-level objects intercept the neutron. 
After the top-level neutron scatterer is identified, the neutron is propagated to 
the front surface of the scatterer if necessary (not necessary if the neutron is 
already inside the scatterer) with appropriate attenuation, and then the ray 
tracing algorithm recurses into itself if the scatterer is a composite. Otherwise 
a point in the forward path of the neutron inside the scatterer will be randomly 
picked, and the neutron is propagated to that point with attenuation. At this 
point a scattering (or an absorption) mechanism of the scatterer is randomly 
picked, and the neutron will be either scattered with its probability adjusted, 
or be absorbed. 

Multiple Scattering 

Multiple-scattering (MS) is naturally supported in MCViNE scattering com¬ 
posites, implemented with a recursive algorithm. For the purpose of this discus¬ 
sion, we differentiate between two types of multiple scattering: single-scatterer 
multiple scattering (SSMS) to describe the multiple scattering within one neu¬ 
tron scatterer and multiple-scatterer multiple scattering (MSMS) to describe 
multiple scattering among neutron scatterers. 


Shape 


Kernels 



Powder-diffraction kernel 
Single -Phonon kernel 
Multi - phonon kernel 


Figure 2: An example of multiple scattering within one scatterer. The incident neutron 
was scattered three times by three different scattering kernels. At each scattering point, the 
original neutron is also propagated out of the scatterer with proper attenuation. Red arrows 
are paths of neutron propagation. Circles highlight the location of scattering. Different 
scattering events are coded using different colors. 


(a) 

(b) 

(c) 

(d) 

(e) 



Figure 3: An example of multiple scattering in a concentric sample assembly. Five (out of 
infinite) possible multiple scattering paths are illustrated. 


Figure 2 shows an example of SSMS in which a neutron gets scattered three 
times before exiting a scatterer. Each time a neutron is scattered inside a homo¬ 
geneous scatterer, the original incident neutron packet is split into two neutron 
packets for computational efficiency. One neutron packet is propagated through 
the scatterer with its probability lowered by attenuation, while the other is scat¬ 
tered by one of the scattering kernels, chosen by a random selection, at a point 
also randomly selected along the forward path of the incident neutron. This 
splitting process repeats for the scattered neutron several times, as illustrated 
in Figure 2, until either the neutron probability is lower than a pre-selected 
limit, or the maximum order of multiple scattering is reached. At this point, 
the neutron that is still inside the scatterer is allowed to propagate out, with 
its probability attenuated appropriately. 

Figure 3 shows an example of multiple scattering in a sample assembly con¬ 
sisting of a “concentric” arrangement of a sample and a sample environment. 
Only five of an infinite number of possible multiple scattering paths are illus¬ 
trated, and the splitting processes are not shown. 
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The multiple scattering algorithm of MCViNE is generic, and it can handle 
more complex sample assemblies, such as non-concentric arrangements, and 
mixing of non-concentric and concentric arrangements. More details about the 
MS algorithm of MCVINE can be found in the supplemental material for this 
paper. 

Sample components in some linear-chain-based MC ray tracing packages can 
support SSMS, but they do not have abstractions similar to composite scatterer 
or scattering kernel. As a result, the SSMS algorithm must be duplicated in 
these sample components, while in MCViNE, implementations of new scattering 
kernels can be added without reimplementing the multiple scattering algorithm. 
McStas supports MSMS partially for a concentric sample assembly by, for ex¬ 
ample, adding a second outer cylinder into the simulation component chain [57]. 
Path (e) of Figure 3 is included this way, but path (c) of Figure 3 is not. Other 
ways of simulating MSMS may be possible with McStas, but these would require 
significantly more effort on the user’s part. 

The multiple scattering algorithm of MCViNE is more comparable to that 
of MSCAT [11] in which the sample and the sample environment are treated 
together in the multiple scattering loop, offering complete treatment of SSMS 
and MSMS. However, MCViNE allows a straightforward increase in complexity 
of the sample and sample environment interactions by using the recursive MS 
algorithm enabled by OOP. MSCAT by default only allows for a few specific 
arrangements. MSCAT and some McStas components contain the ability to 
appropriately transform the angular scattered distribution off of the sample so 
as to only simulate neutrons that will impact the detector, such optimization is 
not yet available in MCViNE. 

Ray tracing in a detector system 

Ray tracing of a neutron through a sophisticated detector system in instru¬ 
ments such as ARCS and SEQUOIA, where flat detector packs are arranged in 
an approximately cylindrical arrangement, illustrates another strength of the 
software design of MCViNE. MCViNE takes advantage of the hierarchical rep¬ 
resentation for neutron scatterers, only in this case the elemental homogeneous 
neutron scatterer is the 3 He detector tube that intercepts neutrons and records 
them. MCViNE reuses the code for ray tracing in a composite scatterer for 
simulating 3 He detector systems, and the new code needed is a scattering kernel 
for the 3 He material that takes into account gas absorption. The ray tracing 
through a cylinder takes care of the parallax effect of the detector tube. When 
a neutron is sent to a detector system shown in Figure 1(c), for example, the 
generic ray tracing algorithm for composite scatterers first checks whether the 
top level composite scatterer is penetrated by the neutron. If so, all constituents 
of the composite scatterer, i.e. the detector packs, are examined to determine 
which of them intercepts the neutron. Unless a neutron traverses a gap be¬ 
tween detector packs, the detector pack is identified and then its constituents, 
the 8 detector tubes, are examined for neutron detection. The path of a neu¬ 
tron through the detector tube is then computed by ray tracing of the neutron 
through a cylinder, and a MC sampling picks a point in the path for the neutron 
to be absorbed. The position in the tube is used to calculate a unique pixel iden- 
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tifier known as the Pixel ID, according to the scheme readable by the Mantid 
software framework [58, 59]. Additionally, the appropriate weighting multiplier 
for neutron probability (computed from absorption probability depending on 
the 3 He pressure and the length of the neutron path through the tube) and the 
time-of-flight are computed for the neutron to be recorded as a detector event in 
a “virtual detector electronics device”. Our hierarchical approach to detectors 
allows for the addition of more details, including details of the charge cloud and 
the wire if deemed necessary. 


3. Examples 

This section presents examples of MCViNE simulations performed for ex¬ 
periments on the ARCS, SEQUOIA, and HYSPEC instruments at SNS. These 
simulations were performed on one of the SNS data analysis clusters. A typical 
analysis cluster has 64 Intel or AMD CPU cores at ^3GHz, and simulations 
usually run in parallel using 10 cores. Every simulation consists of 4 steps: 

Beam simulation. The incident beams on the sample for ARCS, SE¬ 
QUOIA, and HYSPEC instruments were simulated. The MCViNE simulation 
models were derived from the McStas[12] instrument definitions used in the de¬ 
sign phase of the instruments [45, 48, 60, 61]. In this work, 10 9 neutron events 
emitted from the moderator were included in all the simulations, and the neu¬ 
trons at the sample position were saved and reused in the next step. A typical 
parallel run took about one hour. 

Sample scattering. The neutron packets saved in the previous step were 
sent to a SampleAssembly component. Typically 10 8 neutron packets were 
simulated in a run for any powder sample presented in this section, and it took 
minutes to hours to finish, depending on the complexity of the sample assembly. 
The scattered neutrons were then saved. 

Detector interception. Each neutron scattered by the sample was pro¬ 
cessed by a DetectorSystem component and an event was recorded in an event¬ 
mode NeXus Hie if it intercepted a detector tube. In ray-tracing through a 
detector system, the detector tube in which the event was detected was lo¬ 
cated using a hierarchical set of detector tubes as described previously, while 
the exact location and time-of-flight were determined by a MC selection. The 
running time of this step depends on the number of neutrons hitting the detec¬ 
tors, which in turn depends on the number of neutron packets simulated in the 
previous step, as well as the user choice of multiple-scattering. Typically, 10 8 
neutron packets took half an hour. 

Reduction. The NeXus data generated in the previous step were reduced 
using Mantid[58, 59]. The only difference between the simulated and measured 
data is that the intensities in the simulations are computed as the sum of the 
probabilities of all packets arriving in the bin of interest, while those in the data 
are total event counts. The data reduction workflow is therefore identical, and 
uses the same code base for both the simulated and the measured data. This 
step typically took 10s of minutes to finish for a powder measurement. 
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The goal of this four-step simulation process was to reproduce the experiment 
(including data reduction) by a simulation of high fidelity. 

This 4-step simulation workflow produces output files such as the simu¬ 
lated scattered neutrons, the simulated event-mode NeXus Hie, and the reduced 
I(Q 1 E) file. 

The simulation examples here make use of the following scattering kernels 
that are described briefly in the Supplemental Material: 

• incoherent elastic scattering 

• coherent elastic scattering from powder sample 

• incoherent inelastic single-phonon scattering 

• coherent inelastic single-phonon scattering from a powder sample 

• multi-phonon scattering 

• scattering from a dispersion surface where the dispersion relation and 
the dynamical structure factor are described by analytical functions of 
momentum transfer vector Q 

Three examples are presented in the following subsections, two modeling 
vibrational excitations, and one modeling magnetic excitations. The simplest 
sample, a vanadium plate, is presented first in Subsection 3.1, where only in¬ 
coherent scattering kernels are used. Coherent phonon scattering kernels for 
powder samples are introduced for an aluminum sample in Subsection 3.2. Both 
Subsection 3.1 and 3.2 demonstrate the ability to easily turn on and off different 
scattering mechanisms, allowing researchers to gain a better understanding of 
their contributions. Subsection 3.3 presents a simulation of a measurement on a 
single-crystal K 2 V 3 O 8 sample, for which a dispersion-surface scattering kernel 
is used. 

In addition to these three examples, MCViNE simulations of Uranium Ni¬ 
tride (UN) measurements performed on the ARCS and SEQUOIA instruments 
provide an excellent example of the capabilities presented here [9]. In this case, 
the UN sample exhibited particularly strong multiple scattering due to its size. 

MCViNE simulations were able to reproduce well the multiple scattering, iden¬ 
tified weak scattering from accoustic phonons, and showed that the binary solid 
model is a good explanation of the temperature-dependent broadening of the 
modes of the quantum harmonic oscillator [9] . 

3.1. Vanadium 

Here we present experimental and simulated inelastic spectra for a 50 mmx 50 mmx 1.2mm 
vanadium plate sample in the ARCS instrument. A quick calculation using the 
total scattering cross section for V shows that such a sample is a 6 % scatterer. 

Vanadium, being the regular calibration standard in neutron scattering experi¬ 
ments, was chosen as the first example for its simplicity (incoherent scattering 
cross section is much larger than that of coherent scattering, and hence the 
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Figure 4: I(Q,lj) plots of vanadium inelastic spectra obtained from experiment and simula¬ 
tions for a vanadium plate at room temperature in the ARCS instrument, (a) Experiment 
(b) Simulation without a multi-phonon kernel. Multiple scattering was turned off. (c) Simu¬ 
lation with a multi-phonon kernel. Multiple scattering was turned off. (d) Simulation with a 
multi-phonon kernel. Multiple scattering was turned on. 
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Figure 5: Energy spectra integrated over Q range (8,12) inverse angstrom, obtained from 
experiment and simulations for a vanadium plate placed in the ARCS beam. The intensity 
axis is in log scale, (a)-(d) same as Figure 4 


scattering is isotropic). In the experiment, the sample was approximately per¬ 
pendicular to the beam, and the incident energy was tuned to 117 meV using a 
Fermi chopper with 1.5 mm slit spacing and 550 mm radius of curvature, rotat¬ 
ing at 600 Hz. In the simulation, the instrument parameters used matched the 
experimental ones. The simulation sample assembly contains only one homoge¬ 
neous scatterer for the vanadium plate, which was tilted 96.6 degrees from the 
beam direction. Different scattering kernels were used for different simulations, 
but one incoherent elastic kernel and one single-phonon incoherent inelastic 
kernel were included for all. Incoherent scattering cross section of 5.08 barns 
and absorption cross section of 5.08 barns [62] were used. All phonon-releated 
scattering kernels use the phonon density of states (DOS) calculated from a 
Born-von Karman (BvK) model [63], which used force constants originally re¬ 
ported in Ref. [64], and tabulated in the Landolt Bornstein series of material 
properties [65]. Only one universal scale factor was applied to the intensities of 
all simulated spectra to match the experimental data. 

Shown in Figure 4 are I (Q, ui) plots. The top panel is the experimental result. 
Panel (b) was simulated without either multi-phonon scattering or multiple¬ 
scattering. Panel (c) was simulated with multi-phonon scattering but multiple¬ 
scattering was turned off. Comparing (b) and (c), we can see how multi-phonon 
scattering contributes to intensities at high Q (>~ 1(M -1 ), especially notice¬ 
able at high energy transfer above the single-phonon cut-off energy of ~35 meV. 
Panel (d) was done with both multi-phonon and multiple-scattering contribu¬ 
tions. Comparing (c) and (d) shows how the effects of the multiple scattering are 
small but observable, particularly at low Q values (<~ 2A -1 ). As the plate is 
only a 6% scatter, one would expect the multiple scattering to be less than 0.4%, 
consistent with the observation. It also shows that there is less Q-dependence in 
multiple-scattering than in multi-phonon scattering. Figure 5 shows the energy 
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spectra integrated over momentum transfer range of 8 to 12 A -1 in log scale. 
Similarly, we found that multi-phonon scattering contributes more substantially 
than multiple scattering in this Q range. 

To examine the agreement between the simulated and the experimental spec¬ 
tra, constant momentum cuts and constant energy cuts of the experimental data 
(Figure 4(a)) and the simulated data with the most complete collection of scat¬ 
tering kernels (Figure 4(d)) are presented in Figure 6 and Figure 7, respectively. 
The simulated data in general agree very well with the experimental data, across 
the measured dynamic range and the intensity range, showing features such as 
the asymmetric line shape of the elastic peaks, the detailed balance, and the In¬ 
dependent phonon-induced energy-loss profiles in the constant-momentum cuts, 
and the dark-angle dips in the constant-energy cuts. Also in the constant en¬ 
ergy cuts (Figure 7), spikes in intensities (e.g., at E = —20 meV, Q ~ 10 A -1 , 
and at E = 20 meV, Q ~ 9 A -1 ) result from insufficient filling of the displayed 
Q, u! histogram bins by the detector pixels. While additional analysis may be 
able to correct for this effect, it is instructive that the detailed modeling of the 
ARCS detector array almost exactly matches the measured data, while typical 
MC simulations would miss this effect due to simplified detector models. 

The uncertainties of this simulation arise from the following factors: devia¬ 
tion of the simulated incident beam from the actual neutron beam at the ARCS 
instrument; differences between the sample shape specification in simulation 
and the actual sample shape, such as its dimensions, and its orientation relative 
to the incident beam; the uncertainties in the physical and geometrical prop¬ 
erties of the detector packs in the detector system, especially their positions 
and orientations; the deviation of the sample temperature in simulation from 
the actual temperature; and the uncertainties in the sample material such as 
impurities and contamination, neglect of Bragg reflections in simulations, and 
uncertainties in the phonon density of states (DOS) of vanadium. 

The simulation of the incident beam is in part validated by comparison to 
the first beam monitor, and scattering from 2,5-diiodothiophene (C 4 H 2 I 2 S) [47]. 
This partial validation is corroborated by the agreement observed in the elastic 
peaks of I(-E’) spectra (Figure 6 ). The sample dimensions and its position and 
orientation influence the position of the elastic peaks, and the position and the 
width of the dark angle, which all show excellent agreement between experi¬ 
mental and simulated data. Uncertainties in positions and orientations of the 
detector packs can result in incorrect mapping from pixels to Q, E coordinates. 
The fact that the simulated data show the same artifacts at the low-pixel-count 
detector area as the experimental data indicates that these errors are minimal. 
The deviation of experimental sample temperature from the simulation value 
is not expected to be larger than 10 Kelvin, which is consistent with excellent 
match in detailed balance between the simulation and the experiment shown in 
Figure 6 . An examination of Figure 7(b) shows a series of sharp peaks, which 
arises from the Bragg reflections in the sample that were not included in the 
simulations. These may account for some differences in the elastic peaks in Fig¬ 
ure 6 as well. The excellent agreement in all the other cuts of Figure 7 shows 
this does not affect the spectroscopic results. The remaining factor is the un- 
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Figure 6: Constant momentum cuts of spectra as shown in Figure 4 at different Q val¬ 
ues: (a) 8 and (b) 11 A -1 . The experimental cuts were taken from Figure 4(a), while the 
simulated cuts were taken from the best simulation incorporating all scattering kernels and 
multiple-scattering, Figure 4(d). The spectra show decreasing intensities in elastic peaks and 
increasing relative strength of inelastic scattering on both sides of the elastic peaks, with 
larger momentum transfers. The deviations of the simulated data from the experimental data 
(A I = /exp — / s j m , s i m il ar thereafter) shown in the lower panels are small and are within 
error bars. The slightly larger deviations around elastic peaks may be attributed to slight 
mist match in determination of elastic lines, the fact that the coherent scattering is neglected 
in the simulations, as well as uncertainties in moderator profile, and in guide and chopper 
simulations. 
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Figure 7: Constant energy cuts of spectra as shown in Figure 4 at different E values: (a) -20 
(b) 0 (c) 20 and (d) 40 meV. The experimental cuts were taken from Figure 4(a), while the 
simulated cuts were taken from the best simulation incorporating all scattering kernels and 
multiple-scattering, Figure 4(d). The elastic scattering in (b) weakens with higher Q due to 
the Debye-Waller factor; the inelastic scattering in other panels strengthens with higher Q : 
for panels (a) and (c), this trend is dominated by the one-phonon scattering proportional to 
the Q 2 factor in addition to the same Debye-Waller factor, while for panel (d) for which the 
energy transfer of 40 meV is beyond the single-phonon energy cut-off, it is dominated by the 
two-phonon scattering proportional to the Q 4 factor. The dip shown in every panel is typical 
of the plate sample geometry and a result of absorption along the sample width (the dark 
angle). The deviations of the simulated data from the experimental one (shown in the lower 
panels) are in general smaller than error bars. The peaks in the deviation curve at the lower 
panel of (b) are beyond error bars, but they match positions of powder diffraction peaks as 
indicated by the vertical dotted lines, and most likely can be accounted for by the fact that 
the coherent scattering from V are neglected in the simulations. 
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Figure 8: Experimental and best-simulated energy spectra obtained by integrating over Q 
range (8,12) inverse angstrom of the spectra as shown in Figure 4 (a) and (d). The intensity 
scale is expanded to highlight the inelastic spectra. 


certainties in the phonon DOS, which is the main source of uncertainties in this 
simulation. 

The discrepancy between the simulated data and the experimental data, 
relevant to researchers interested in phonons, is easier to observe from Figure 8. 
The double peaks near 20 and 28 meV both appear to shift to the right in 
the experimental data, compared to the simulated one. This discrepancy has its 
origin in the phonon DOS used in the simulation, which is computed from a BvK 
model fitted to an X-ray determination of phonon dispersions [64], Consistent 
with the discrepancies observed here, this BvK model slightly underestimated 
phonon energies when compared to a phonon DOS spectrum measured using 
inelastic neutron scattering (see Figure 10 of [64]). Actually, the uncertainty 
of measured vanadium phonon DOS reported over the years is larger than the 
discrepancy observed here (see, for example, Figure 2 of [66]). 

3.2. Aluminum 

In this example, we present experimental and simulated inelastic spectra 
from a 60 mmx60 mmx4 mm polycrystalline aluminum (1100 alloy) plate in the 
ARCS instrument. A quick calculation using the total scattering cross section 
for A1 shows that such a sample is a 5% scatterer. In the experiment, the sample 
was placed approximately at 135 degrees from the beam, and the incident energy 
was tuned to 80.5 meV using the same Fermi chopper slit package as used in 
the previous example for vanadium, spinning at 480 Hz. In the simulation, the 
instrument parameters used matched the experimental ones. The simulation 
sample assembly contains only one homogeneous scatterer for the aluminum 
plate. Different simulations made use of different combinations of scattering 
kernels. Tabulated A1 cross sections [62] were used. All phonon-related scatter¬ 
ing kernels use phonon energies and polarization vectors computed on a regular 
grid in a Brillouin zone from a BvK model [67]. The broadening of the phonon 
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Figure 9: I(Q,cu) plots of aluminum inelastic spectra obtained from experiment (a) and 
simulations (b)-(e) of an aluminum plate at room temperature in the ARCS instrument. 
Scattering kernels included in the simulations are (b) incoherent elastic and incoherent inelastic 
single-phonon scattering (see note on intensity scaling in text), (c) coherent elastic (powder 
diffraction) and coherent inelastic single-phonon scattering, (d) all kernels included in (b) and 
(c) plus a multi-phonon kernel using an incoherent approximation and (e) all kernels included 
in (d) plus multiple scattering. 
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Figure 10: Energy spectra integrated over Q range (6,10) inverse angstrom, obtained from 
experiment and simulations for an aluminum plate placed in the ARCS beam. The intensity 
axis is in log scale, (a)-(e) same as Figure 9 


modes was not included. Only one universal scale factor was applied to all 
simulated spectra to match the experimental data. 

Figure 9 shows I(Q,to) plots for the aluminum plate. The experimental 
result is given in panel (a) and panels (b)-(e) show simulated data. In (b) only 
the incoherent elastic and incoherent single phonon scattering are included. In 
this plot the maximum intensity for the color scale is reduced by the ratio of 
incoherent/coherent cross sections of aluminum, so one can see the details in the 
plot. The numerical values on the colorbar give an indication of the scaling. In 
(c) only the coherent elastic (powder diffraction) and the coherent single-phonon 
inelastic scattering are included. In (d), all of the kernels in (b) and (c) with the 
addition of a multi-phonon kernel using the incoherent approximation. In (e), 
all of the kernels in (d) are used with multiple scattering turned on. Overall the 
features shown in the experimental data (a) and the simulated data (e) agree 
very well. Comparison of (b) and (c) shows that coherent scattering gives rise 
to more features such as diffraction peaks and phonon dispersion curves. It 
is evident from comparing (c) and (d) that multiphonon scattering increases in 
intensity at higher Q. The most obvious difference in (d) and (e) is in the elastic 
line, which shows that multiple scattering of coherent elastic scattering seems 
to contribute similarly to incoherent scattering in the elastic line. The elastic 
lines in (a) and (e) seem to show that the sample used in the experiment may 
contain traces of an additional phase, most likely from a surface layer of AI 2 O 3 . 
Figure 10 shows energy spectra integrated over a range of momentum transfer 
from 6 to 10 A -1 in log scale. Similar to the vanadium measurement, we found 
that multi-phonon scattering plays a significant role in the observed intensity, 
whereas multiple scattering plays little role in this Q range. The is expected as 
the sample is only a 5% scatterer. 

Similar to the vanadium study in the previous section, we further examine 
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Figure 11: Constant momentum cuts of spectra as shown in Figure 9 at different Q values: 
(a) 3 (b) 5 (c) 7 and (d) 9 A -1 . The experimental cuts were taken from Figure 9(a), while the 
simulated cuts were taken from the best simulation incorporating all scattering kernels and 
multiple-scattering, Figure 9(e). The deviations of the simulated data from the experimental 
data are plotted in the lower panels. The main contribution to the larger deviations around 
the elastic peaks most likely comes from the texture in the sample. 
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Figure 12: Constant energy cuts of spectra as shown in Figure 9 at different E values: (a) -20 
(b) 0 (c) 20 and (d) 40 meV. The experimental cuts were taken from Figure 9(a), while the 
simulated cuts were taken from the best simulation incorporating all scattering kernels and 
multiple-scattering, Figure 9(e). The deviations of the simulated data from the experimental 
one (shown in the lower panels) in general are smaller than the error bars. The discrepancies 
in (b) are from differences in the measured and calculated Bragg peak intensities. These arise 
from the fact that the sample is not a perfect powder. 
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the agreement between the simulated and the experimental spectra by taking 
constant momentum cuts and constant energy cuts of the experimental data 
(Figure 9(a)) and the simulated data with the most complete collection of scat¬ 
tering kernels (Figure 9(e)). These cuts are presented in Figure 11 and Fig¬ 
ure 12 , and they again display strong agreement between the experimental and 
simulated data. For Figure 11, the maximum deviations occur near the elastic 
peak. Similarly, the biggest differences in Figure 12 are observed for the elas¬ 
tic scattering in (b). These indicate that the measured Bragg peak intensities 
vary from the model: an isotropic powder average for Al. This could arise from 
texture, preferred grain orientation, or incomplete sampling of the powder[ 68 ]. 

The uncertainties of this simulation, like the case of vanadium, arise mainly 
from the uncertainties in the phonon data, which is the major input to all 
phonon-related kernels used in the simulation, only in this case the phonon 
data is not just a simple DOS curve. This phonon data input, existing in 
the form of energies and polarization vectors of all phonon branches, was com¬ 
puted for all points on a regular grid inside the first Brillouin zone, from a 
BvK model fit to phonon dispersions measured along the [100], [110], and [111] 
directions [67]. Small changes in the force constants in the BvK model can re¬ 
sult in large fluctuations of the computed phonon energies and polarizations. 
Hence, the discrepancy between the experimental data (Figure 10(a)) and the 
simulated data (with all kernels, Figure 10(e)) was not unexpected. Another 
possible contributing factor is the fact that the broadening of Al phonon modes 
[69, 70] was not taken into account in our simulation. 

3.3. K 2 V 3 Og single crystal 

We now present experimental and simulated inelastic single-crystal spectra 
for a K 2 V 3 O 8 sample [71] measured at the HYSPEC instrument [51], a hybrid 
spectrometer at SNS with a focusing monochromator and a movable detector 
vessel. The sample consisted of 5 co-aligned crystals with an overall cylindrical 
shape approximately 3.8 cm in diameter and 2.5 cm in height. The sample was 
oriented so that at zero degrees of rotation angle its [ 100 ] direction was along 
the beam and its [001] direction pointed upward vertically. Measurements were 
performed at 1.5 K with an incident energy of ~7 meV with a Fermi chopper 
frequency of 180 Hz. The detector vessel was oriented so that neutrons scattered 
from the sample in the horizontal plane were measured in the scattered angle 
range of -75 to -15 degrees. To span reciprocal space, the goniometer angle was 
swept from 40 to -50 degrees in 0.5 degree steps and then from -50 to -28 degrees 
in 1 degree steps. 

The simulation was done for a scan matching the experimental setup. The 
sample assembly contains one homogeneous scatterer for the K 2 V 3 O 8 sample 
with a shape matching the overall shape of the sample in the experiment. Two 
scattering kernels were included in the simulation: one incoherent elastic scatter¬ 
ing kernel to approximate the elastic line, and one dispersion-surface scattering 
kernel for simulating the scattering from the spin-wave. 

In K 2 V 3 O 8 the long-wavelength spin-wave dispersion is well described by a 
result for the quantum (S=l/2) square lattice Heisenberg antiferromagnet: 
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Figure 13: Slices along the Q=[/ilO] axis and the energy axis from (a) experiment and (b) 
simulation, taken in the range Z=(-0.3, 0.3) and fc=(0.93,1.07). 


E( Q) = 2jyjl - 1 q // (1) 

with the dynamic structure factor 

S(Q,E) = ySy((ii-Ii(Q)) ( 2 ) 

V 1 - ^ 

where 7 q // = cos(hn) cos(kn), 2 J = 2.563 [71]. The differential cross section is 

{Jidkf) ~ % |F(Q)|2 [n{E) + 1](1 + cos2 </>)5(Qj E) (3) 

where ki and kf are the magnitudes of the incident and final neutron wave 
vectors, F(Q) is the magnetic form factor for V 4+ [72], n(E ) + 1 is the Bose 
occupation factor, and 1 + cos 2 <f> is a polarization term is the angle between 
Q and the easy c-axis). 

In the treatment of both experimental and simulated data, Mantid was used 
to reduce the measured NeXus file to NXSPE format, a specialized intermediate 
inelastic neutron data format, at each goniometer angle, and then a Python 
program was used to project data in NXSPE files to the four dimensional Q, E 
space. Slices along desired Q directions could then be taken. 
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Figure 14: Constant energy cuts of the slices as shown in Figure 13 in the energy range of 
[1.2, 1.3] meV. A background term and a scaling factor were the only parameters used to fit 
the simulated data to experimental data. The goodness of fit, reduced \ 2 , was 2.56. 


Shown in Figure 13 are experimental and simulated slices along axes of 
Q = [hlO] and energy. The spin wave dispersion in the energy range between 
~0.3 and ~2.1 meV shows good agreement between experimental and simulated 
data. In the experimental data, the noise in the energy range of [0.25, 0.7] is 
higher than the noise at higher energies. This is due to elastic scattering from 
the cryostat, which was subtracted from the data. Both plots show a distinctive 
feature - for the dispersion near [010] the right branch is sharper than the 
left branch, while the dispersion branches near [210] are largely symmetric in 
their broadening. This effect is more clear in Figure 14, where constant-energy 
cuts over [1.2, 1.3] meV are plotted against Q = [hl0] for the experimental 
(circles with error bars) and the simulated (solid line) intensities. Both the 
experimental data and simulated data show that the peak near h = 0.15 is 
much sharper than that near h = —0.15. Such an effect is expected due to 
diffrent focusing conditions of the resolution ellipsoid [73, 74]. The spin-wave 
dispersion above ~2.1 meV (h ~ 0.33, for example) has weaker intensity in 
the experimental data than expected from the simulated data, which could be 
attributed to mode splitting near the zone boundary [71]. 

4. Conclusions 

The MCViNE software package for Monte Carlo neutron ray tracing sim¬ 
ulations is based on modern object-oriented software design that decomposes 
neutron scatterers into a hierarchy of composite and elemental scatterers. This 
software design elegantly solves algorithmic problems including multiple scat¬ 
tering in a sample assembly, and ray tracing in a sophisticated detector sys¬ 
tem. Multiple scattering and instrument resolution are included naturally in 
MCViNE simulations, and simulation of ray-tracing in a complex detector sys¬ 
tem such as those of the ARCS and SEQUOIA instruments is straightforward. 
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Adding, removing, and modifying scattering kernels in the simulations can help 
show the contributions of different types of scattering events to the experimen¬ 
tal data. The examples of scattering from vibrational and magnetic excitations 
presented here and elsewhere [9] demonstrate that such simulations can improve 
our understanding of the underlying physics of neutron spectra. 

Limitations of the MCViNE software package include (1) the complexity 
of sample assembly or detector system can be limited by available computing 
resources (memory and computation time); (2) it currently only supports ho¬ 
mogeneous scatterers, and this means materials with non-uniform distribution 
of defects, for example, can not be simulated without approximations; (3) it 
has a limited but expandable libary of scattering kernels. For example, for 
magnetic systems, it currently only supports spin-wave dispersions that can be 
expressed as analytical functions; (4) MCViNE currently relies on components 
adapted from other MC ray tracing packages to handle many neutron optical 
components such as guides and choppers, and thereby inherits their limitations. 

MCViNE is freely available for the Linux platform. More details about the 
conditions of use and license can be found at http: / /danse.us/trac/MCViNE/wiki/license. 

Details on build and installation, usage, source repository, and user support of 
MCViNE are available in the documentation [32] . Feedback to the MCViNE de¬ 
velopers can be provided through the MCViNE user mailing list http://groups.google.com/group/mcvine- 
users. 
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1. Why Object Oriented? 

Imperative programming is the predominant programming paradigm em¬ 
ployed in MC neutron ray tracing programs. In imperative programming, differ¬ 
ent aspects of programming data and logic are tightly coupled. Object-oriented 
programming (OOP) is less used in scientific programming. Fundamental OOP 
concepts include encapsulation (packing of data and associated functions), poly¬ 
morphism (processing of objects according to type), dynamic dispatch (dynami¬ 
cally deciding which implementation to choose for an operation), and delegation 
(allowing a different object or method to do the work). These methods promote 
abstraction and decoupling of data and programming logic, generally result¬ 
ing in codes that are more extensible, more maintainable, and easier to debug 
[????]. It should be noted that runtime efficiency of a program is not 
the main focus of the OOP paradigm: to achieve polymorphism, for example, 
additional machine codes are generated by compilers to facilitate the dynamic 
dispatch by using virtual tables. OOP codes are also harder to optimize by 
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Figure 1: The MCViNE architecture is divided into three categories for serving different neu¬ 
tron ray-tracing components: simple NumPy-based components, McStas-based components, 
and components built with OOP principles for sample assemblies and detector systems. These 
categories are served by the vertical integration of three layers comprised of Python compo¬ 
nents at the top, followed by python libraries, then CH—h libraries. The foundation of the 
MCViNE architecture is a C++ library with basic neutron objects and math functions. 

the compiler. It is not uncommon for OOP codes to run slower than programs 
written with imperative programming. However, this inefficiency can be over¬ 
come by tremendous improvements in both developer productivity and software 
sustainability. 

2. Layered Architecture of MCViNE 

Figure 1 depicts the architecture of the MCViNE package. Two layers of 
libraries exist in C++. The bottom layer is a basic C++ library that defines 
neutrons and some basic math facilities such as 3D vectors. Above this layer 
are three different types of C++ libraries: 1) A NumPy[? ] array adapter 
that allows a group of neutrons to be manipulated easily as a NumPy array 
at the Python level. 2) A McStas adapter C library that provides support 
for MCViNE-wrapped McStas components. 3) An advanced C++ library that 
supports the concept of neutron scattering composite and scattering kernels, a 
unique feature of MCViNE. Above the three C/C++ libraries, a layer of Python 
helper classes and methods are implemented to provide a Python interface to 
the underlying C++ functionalities. 
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Three types of Python components sit on top of the Python adaptation layer, 
including: 1) NurnPy based neutron components in which a group of neutron 
packets can be manipulated as a NumPy array. This makes it extremely easy to 
create simple neutron components, and to create prototypes of more sophisti¬ 
cated neutron components. 2) MCViNE-wrapped, shared-library-based McStas 
components in which a MCViNE facility automatically takes McStas compo¬ 
nents, compiles them into shared libraries, and binds them to Python. Unlike 
simulations in McStas, MCViNE components are not compiled to a monolithic 
executable. 3) Generic scatterer components that support composite scatterers 
and scattering kernels, which are described in detail in the main text. 

3. Multiple Scattering Algorithm 

To implement the multiple scattering algorithm in MCViNE, the Compos- 
iteScatterer has a method “scatterM”. It takes one neutron as input, and outputs 
all scattered neutrons, which includes the original neutron with its probability 
lowered by attenuation. The pseudo-code (in a python-like syntax) of scatterM 
is 

def scatterM(scatterer, neutron): 

# initialize array of neutrons to be scattered 

toscatter = [neutron] 

nms = 0 # multiple scattering order 

out = [] # initialize array of scattered neutrons 

# main loop 

while toscatter is not empty and nms < scatterer .max.ms : 

# init array of scattered neutrons 

scattered = [] 

# inner loop over all neutrons that need scattering 

for neutron in toscatter: 

# if neutron probability is low, skip 

if neutron.probability < scatterer.min_prob: 
continue 

# if the neutron is leaving the scatterer , save it 

if neutron is exiting: 
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out.append(neutron) 
continue 

# call interactM_pathl to compute all neutrons scattered 

# from the first encouter of the neutron with the scatterer 

scattered += scatterer.interactM_pathl(neutron) 

continue 

# scattered neutrons may need to be scattered again 

toscatter = scattered 

# increment multiple scattering order 

nms += 1 

continue 

# propagate remaining neutrons 

if toscatter is not empty: 

out += propagat_out_with_attenuation(toscatter , scatterer) 

# return scattered neutrons 

return out 

This code contains two loops. The outer loop (the while loop) keeps track of two 
variables, “toscatter” (the neutrons to be scattered), and “nms” (the multiple 
scattering order). For every neutron in the “toscatter” variable whose forward 
path intersects the scatterer, the inner loop computes a set of scattered neutrons 
by calling the method interactM_pathl. This set of scattered neutrons (“pathl- 
scattered-neutrons” for future reference) may be scattered zero, one, or more 
times by the scatterer, but each of them has a (multiple-scattering) path that 
is always inside the first constituent scatterer that it encounters (“scattererl”). 
There can be only one exception: if the original position of the neutron is 
outside the scatterer, the path segment that leads the neutron from its original 
position to its entry point of scattererl can be in vacuum or air. The computed 
neutrons are accumulated into the variable “scattered” which later is renamed 
to “toscatter” for the next round of this inner loop. The inner loop skips over 
neutrons that have probabilities lower than a threshold selected by the user. The 
outer loop keeps running the inner loop until either all the neutrons scattered 
in the inner loop leave the scatterer (no more interceptions), or the maximum 
multiple scattering order, as set by the user, is reached. 
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The “interactM.pathl” method of a composite scatterer takes one neutron 
as input and computes “pathl-scattered-neutrons”. Its implementation finds 
the first scatterer “scattererl” that intercepts the neutron, transforms the neu¬ 
tron to the coordinate system of scattererl, calls the “interactlVLpathl” method 
of scattererl, and finally transforms the scattered neutrons back into the coor¬ 
dinate system of the host scatterer. Scattererl could be a composite scatterer 
itself, in which case the call of the “interactlVLpathl” method would recurse 
into itself. If scattererl is an elemental scatterer, the elemental scatterer’s inter¬ 
actlVLpathl method will be called, hence stopping the recursion. The selection 
of the interactM_pathl method is dynamically done by C++, thanks to the 
polymorphism. 

The only elemental scatterer implemented so far is the HomogeneousScat- 
terer. Its implementation includes an interactM_pathl method that can be 
illustrated by the following pseudo code: 


def interactM_pathl(hs, neutron): 

# input s : 

# hs: homogeneous scatterer 

# neutron: input neutron 

# output: 

# scattered neutrons . these neutrons are scattered zero , one , or 

# more times by "hs", but the paths of those neutrons are always 

# inside "hs" except the leading path to the entry point. 

# 

# initialize the neutron to be scattered 

toscatter = neutron 

nms = 0 # multiple scattering order 

out = [] # initialize array of scattered neutrons 

# main loop 

while toscatter is not null and nms < hs . max_ms : 

# if neutron probability is low, skip 

if toscatter.probability < scatterer.min.prob: 
toscatter = null 
break 

# let neutron interact with the scatterer once 
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hs.interactMl(toscatter) 


transmitted, scattered = 

# transmitted neutron should be saved 

out.append(transmitted) 

# if the neutron is not inside this scatterer 

# done . 

if not isinside(scattered , hs ) : 
out.append(scattered) 
toscatter = null 
break 

toscatter = scattered 

# increment multiple scattering order 

nms += 1 
continue 

# remaining neutrons 

if toscatter is not null: 
out.append(toscatter) 

# return scattered neutrons 

return out 


The pseudo-code is mostly self-explanatory, and the only thing that deserves 
explanation is the method interactMl of a HomogeneousScatterer. This method 
takes one neutron as input, and computes two output neutrons, transmitted and 
scattered. The transmitted neutron is the forward-propagated neutron with its 
probability lowered by attenuation along its first segment of its forward path 
through the scatterer. The scattered neutron is computed by randomly picking 
one point along the first segment of the forward path of the neutron, propagating 
the neutron to that point with attenuation, selecting a kernel and calling the 
kernel’s scatter method to change the state of the neutron, and finally adjusting 
the probability of the neutron. 

4. Scattering kernels used in the examples 

Incoherent elastic scattering is simply given by [? ] 



da Ui 
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where exp(— 2W) is the Debye-Waller factor, and (Ji nc is the incoherent cross 
section. 

Coherent elastic scattering from powder sample. The total cross 
section for scattering neutrons into a Debye-Scherrer cone is [? ] 


°v l-M r/ )l 2 

Vo 4 sin | , 

2 Tf—T 

The sum is over reciprocal lattice vectors of the same length. A is the wavelength 
of the neutron, 6 is the scattering angle, and Fn is the nuclear structure factor. 

Incoherent inelastic single-phonon scattering is given by [? ] 


- = bV ,z(£) 

dfldEf J ±1 4tt h 2 M v ; E 


d 2 a \ 


coth 


/ hco \ 
\2k B T) 


± 1/2 


( 2 ) 


where ki and kf are the momentum of the incident and outgoing neutrons, re¬ 
spectively, Q is the momentum transfer, E is the energy transfer, jcoth ( ^ ± 1 j /2 

is the thermal factor, and Z(E) is the phonon density of states as a function of 
phonon energy E, which is the main input of the kernel. 

Coherent inelastic single-phonon scattering from a powder sample 

is given by [? ] 


(*) 

\dQ,dEf ) 


mc±l 


cr C oh kf (27r) 3 
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where e s is the polarization of the phonon mode, E s is the energy of the phonon 
mode, and r is a reciprocal lattice vector. The total cross section for a phonon 
mode of energy E at Q can be deduced as 
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The main input for this kernel is the energies and polarization vectors of phonon 
modes in a Brillouin zone. 

Multi-phonon scattering is given by [? ] 
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The main input is the phonon DOS. 

Scattering from a dispersion surface E = E( Q) with a scattering 
function of S = S( Q) is given by 
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The main inputs are analytical functions E{ Q) and S( Q) as strings. 
Here is an example of an xml specification for this type of kernel: 

<E_vQ_Kernel 

E_Q = "2.5* sqrt(1-(cos(4.4*Qx)*cos(4.4*Qz))~2) " 

S_Q="1" 

/> 







